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Abstract 



^: 

^ , This paper addresses the behavior of large systems of heterogeneous, globahy coupled oscillators 

"^ ' each of which is described by the generic Landau-Stuart equation, which incorporates both phase 

and amplitude dynamics of individual oscillators. One goal of our paper is to investigate the effect 
of a spread in the amplitude growth parameter of the oscillators and of the effect of a homogeneous 
P^ I nonlinear frequency shift. Both of these effects are of potential relevance to recently reported 

experiments. Our second goal is to gain further understanding of the macroscopic system dynamics 
at large coupling strength, and its dependence on the nonlinear frequency shift parameter. It is 
proven that at large coupling strength, if the nonlinear frequency shift parameter is below a certain 



> 

Q>^ ' value, then there is a unique attractor for which the oscillators all clump at a single amplitude 



and uniformly rotating phase (we call this a single-cluster "locked state"). Using a combination of 
analytical and numerical methods, we show that at higher values of the nonlinear frequency shift 
parameter, the single-cluster locked state attractor continues to exist, but other types of coexisting 
attractors emerge. These include two-cluster locked states, periodic orbits, chaotic orbits, and 



C^ . quasiperiodic orbits. 



PACS numbers: 



Systems of coupled oscillators occur in a very wide variety of applications. 
Often interaction between the dynamical evolution of the oscillator phases and 
amplitudes is an important issue. The simplest model of such dynamics is that 
of a globally coupled system of Landau-Stuart equations |ll[|. While this system 
is very basic, due to the large space of possibilities for its parameters and their 
probability distribution functions, a complete understanding of the system is 
lacking. Here, motivated by recent experiments [2|, we consider parameter de- 
pendences not previously investigated. We also investigate the reason for the 
common occurrence of "locked states" (constant amplitude and sinusoidal oscil- 
lation) observed in previous studies when the coupling is large, why non-locked- 
state attractors occur at sufficiently large values of the nonlinear frequency shift 
parameter, and what types of attractors can occur at large coupling and large 
nonlinear frequency shift. 

I. INTRODUCTION 

The interaction of many coupled dynamical units is a common theme across a broad range 
of scientific disciplines. Within this general theme, the issue of determining conditions for 
the emergence of macroscopic cooperative behavior and of determining the nature of this 



behavior is of centra! 
junction circuits 71, 



importance [jj, |3|] . Examples include coupled lasers 
8|, interacting yeast cells 



-E 



J]-[6|, Josephson 



ipl, pacemaker cells in the heart 11|. 



pedestrian induced oscillation of footbndge. mM, chenJcally .eactn.g sy.ten. U 
circadian rhythms [16|], and many others. 

A very useful simplified framework for beginning to understanding phenomena observed 
in these situations is the phase oscillator description. In the phase oscillator description 
the dynamical units are assumed to be oscillatory with fixed amplitude. Thus, the state 
of each individual unit can be specified solely by a phase angle 6, and the evolution of 
oscillator i is taken to be determined by its present state 6i and by the states 6j of the 
other oscillators (j 7^ i). The simplest model of this type was originally put forward by 
Kuramoto in 1975 and has proven to be an extremely useful paradigm for understanding 

n n 

this general type of system [17|-[21[. In addition, the Kuramoto model has also served 
as a basis for formulating related phase oscillator models appropriate to a wide variety of 



situations (e.g., see Ref. |22|). The basic question addressed by the Kuramoto model is 



the competition between the synchronizing effect of coupling and the desynchronizing effect 



22| 



of different natural frequencies of the individual oscillators. The principal result [17 1 
coming from the solution of the Kuramoto model is that, in the limit of a large number of 
oscillators, this competition is resolved by a transition, whereby there is a critical coupling 
strength below which the oscillations of individual oscillators occur with random phase and 
there is no macroscopic population-wide oscillation, while above which oscillators begin to 
develop phase coherence, and globally-averaged population-wide oscillation emerges. 

A drawback of the phase oscillator description is that, by its definition, it excludes the 
effect of amplitude dynamics and the possible coupling of amplitude dynamics with phase 
dynamics. Another useful oscillator model incorporating both amplitude and phase dynam- 
ics is based on the normal form of an isolated system near a Hopf bifurcation, 

dz 

— = {a + iuj)z-{(5 + i-i)\z\^z, (1) 



also referred to as a Landau-Stuart oscillator 



a 



In ([T]) 2; is complex with | z \ and the angle of 



z representing the amplitude and phase of the oscillator. The real parameter a is the linear 
amplitude growth rate of the oscillations, with a > for growth (and a < for damping). 
The Hopf bifurcation occurs as a passes through zero. The other real parameters a;,/3,7 
respectively characterize the small-amplitude natural frequency of the oscillator, and the 
finite amplitude nonlinear shifts of the small amplitude growth rate and frequency. The 
bifurcation is supercritical if /3 > (the nonlinear term saturates growth) and subcritical 
(hysteretic) if /3 < (the nonlinear term enhances growth). Here we will only deal with the 
supercritical case [in the subcritical case, if a > 0, orbits typically go far from 2; = 0, thus 
invalidating the expansion resulting in ([1])]. For a < 0, Eq. ([1]) has as its stable solution 
z = 0. For a > 0, z = is unstable, and ([1]) results in an attracting limit cycle attractor. 



ja 
exp 



w - ^ 1 t + ^n 



^,^.^„. (2) 

which traces a circular orbit about the origin of the complex z-plane. In general, the normal 
form oscillator parameters (a,a;,/3,7) derived from the physical system under study will 
depend on some set, p = {j)'^^\p^'^\ ■ ■ ■ ,p'^"))-^, of physical parameters for that system. That 
is, [cu,a,/3,7] = [w(p),a(p),/3(p),7(p)]. 



We are interested in the situation, also studied in Refs. [23|]-[29|, where many oscillators 
of the form of Eq. ([1]) are coupled together and where each such oscillator (indexed by a 
subscript i = 1,2,--- ,A^ ^ 1) may have a different parameter set. That is, if oscillator i 
has parameter set pj, then 

[uji,ai,/3i,'ji] = [w(pi),a(pi),/3(pi),7(p,)]. (3) 

If the value of p is regarded as assigned randomly from oscillator to oscillator according to 
some probability distribution function (pdf), then that will induce a corresponding pdf G 
of the parameters [u,a, /3,'j], such that 

G{uj,a, (3,'y)dujdad(3d'y (4) 

represents the fraction of oscillators with parameters u, a, (3, 7 in the range u E [u,ijJ + du], 
a E [a,a + da], /3 G [f3, (3 + df3], 7 G [7, 7 + d'j], and apphcable in the limit A^ — )■ 00, where 
A^ is the number of oscillators. 

Considering this general problem, one would like to know how the system behavior de- 
pends on the distribution function G{u!, a, (3, 7). However, as G is a distribution in the four 
variables a;,a,/3,7, this is clearly too big a problem to address in full generality. Here we 
will pursue a more modest program. In particular, the questions we address are partly mo- 
tivated by the experimental work in Ref. pi]: (i) what is the effect of spread in a allowing 
the simultaneous presence of dead (a < 0) and active (a > 0) oscillators in the uncoupled 
state, and (ii) what is the effect of a nonlinear frequency shift 7 (for simplicity, we treat the 
oscillators as all having the same 7), and iii) how we can understand certain types of sim- 
ple nonlinear behavior often observed in these systems when the coupling strength between 
oscillators is large? 

II. FORMULATION, BACKGROUND AND OUTLINE 

We assume /3 and 7 are the same for all oscillators, Pi = /3 and 7^ = 7. Furthermore, we 
scale /3 to one by a proper normalization of Zi {zi — )■ Zi/^//3). Thus 

Gioj, a, (3, 7) = Giu, a)5{P - 1)5(7 - 7)- (5) 



If Lo and a are uncorrelated in their variation from oscillator to oscillator, then G is of the 
form 



G{uj,o.) = g{uj)h{a). 



(6) 



In what follows we assume that Eq. (I6l) holds 30|, and that g{uj) is symmetric and mono- 
tonically decreasing with respect to its maximum value, which we can take to be located at 
w = (if the maximum of g{uj) occurred at some non-zero value, w = a), then the location 
of the maximum can be shifted to zero by the change of variables z = z'e*"^*, u' = u — u). 

For ([T]) with /3j = 1, 7, = 7 and (jS]) specifying our ensemble of uncoupled oscillators, we 
now proceed to globally couple these ensemble members through a mean field, (z), 



— ^ = (oi + ioJi)zi - (1 + i'y)\zi\'^Zi + V{z), 

N 



(^)-^E 






(7a) 
(7b) 



where the parameter F measures the strength of the coupling and is assumed real and pos- 
itive, F > (some previous studies have considered complex coupling constants, e.g., Refs. 
371]- 40|). We will sometimes refer to {z) as the "order parameter" because whether or not 
there is global collective behavior for A^ — )■ 00 corresponds to whether |(2;)| > or [z) = 0. 
See Refs. |23|-[29| for previous related work on large coupled systems of Landau-Stuart equa- 
tions. In many of these previous works 23|- 27|], the coupling term is written as K{{z) — Zi) 



in place of T{z). This choice is simply related to ours by the transformation «« = dj — K, 
r = K. We prefer the parametrization in Eq. ([7]) because one of our principal motivations 
will be experiments [2(] where it can be plausibly argued that quantities analogous to F and 
the average value of a^ (denoted a) can be varied essentially independently. More gener- 
ally, in real large coupled oscillator systems familiar to us, coupling between the oscillators 
typically results from physical processes distinct from those determining the properties of 
the individual oscillators (as in Ref. |2|), and the parametrization in Eq. ([71) is therefore 
the appropriate one. Use of the form ([7]) (rather than ([8]) below) will be important for our 
considerations of the large coupling limit in Sec. IIX[ In addition, in Refs. 23|- 27| it was 
considered that dj was the same positive constant for all i, &{ = a, and furthermore that 



7 = 0. Parameter and time normalizations were then chosen to transform d to 1, yielding, 
in place of ([7]) 

dz- 

-^ = {iuji + l-\zi\'')z, + K{{z)-Zi), (8) 

where the coupling parametrization form K{{z) — Zi) was used. We, however, will be inter- 
ested in the effect of a spread in a^ with the possibility of the simultaneous occurrence of 
positive and negative ctj for different i, and also in the effects of nonlinear frequency shift 
7 7^0. 

One motivation for this study is the recent paper, Ref. |2|], which describes experiments 
in which many (~ 10^/cm^) specially designed small porous particles are continuously and 
rapidly mixed in a catalyst-free Belousov-Zhabotinsky reaction mixture. The catalyst for 
the reaction is immobilized on the small porous particles, each of which can potentially serve 
as an effective chemical oscillator. Oscillations in the chemical states of the particles are 
visualized as the color of the particles oscillates between red and blue. The particle density 
serves as a parameter analogous to our coupling constant F, while regulation of the stirring 
rate effectively provides a control analogous to control of the mean oscillator growth rate, 

a = ah{a)da. (9) 



Because the process by which the particles are prepared is not perfect, it is expected that 
there will be substantial spread in their parameters, and in particular in u and a. These 
spreads are of particular interest because: (i) spread of oscillator frequencies is the essential 
feature leading to the transition from incoherently oscillating units to macroscopic oscillation 
in the Kuramoto model, and (ii) the parameter a determines whether individual particles, 
when uncoupled, oscillate (a > 0) or do not oscillate (a < 0). In the case a < the attractor 
for Eq. ([1]) is the fixed point z = 0, often referred to as "oscillator death". With reference to 
point (ii), because of the spread in a, in some range of stirring rates, we can expect a situation 
like that shown schematically in Fig. [1], which depicts an uncoupled oscillator growth rate 
pdf h{a) yielding substantial fractions of the particles in the oscillating and dead states. 
As a increases from very negative values, (— «) ^ Sa) (analogous to low stirring rates in 
the experiment), to very positive values, a ^ 6a (analogous to high stirring rates), there is 
a continuous transition from predominantly dead to predominantly oscillatory dynamics of 



the uncoupled oscillators. Another notable feature of these experiments is that the collective 
coherent frequency of oscillation exhibits a marked dependence on the oscillation amplitude 
through its dependence on the density of the porous particles at fixed stirring rate (e.g., the 
third panels in Fig. 2(a) and 2(b) of Ref. J2l]). This is a strong indication that the nonlinear 
frequency shift 7 plays a significant role. It is notable that Ref. [21] developed a set of chemical 
rate equations that, when solved numerically, yield good agreement with the experiments. 
While this is a singular achievement, we are interested in obtaining additional understanding 
of the processes involved and in determining if it is generic. To the extent that qualitative 
behavior of our Landau-Stuart model mimics behavior observed in the particular experiment 
in Ref. |2| , the typicality of the observed phenomena is strongly implied. Furthermore, if the 
above agreement holds, then any analytical results obtained for the Landau-Stuart model 
may lead to further understanding of these experimental phenomena. Thus it is our desire 
to employ the generic coupled Landau-Stuart model, Eqs. ([7]), to explore and understand 
the nature of the interplay between frequency spread, growth rate spread and nonlinear 
frequency shift. In this connection, it is worth noting that our work may be applicable to 
other experiments. Indeed, as described in Ref. 21, the chemical experiment was, at least 

nn 

partly, intended to mimic observed oscillator quorum-sensing in yeast populations [9|, |10J. 
In addition, the basic stability analysis technique used here (Sect. IIIip is similar to that 
originally introduced in Refs. 26| and 27| can also be applied to other amplitude / phase 
oscillator systems, such as the laser system considered in Ref. pj. 




FIG. 1: Schematics of h{a) 



We now giv e a brief review of the most important papers 23|]- 29| related to our work. 



References 



23|-[27| considered Eq. ([8]) (all oscillators have identical aj and 7j = 0) and 



examined the behavior as a function of the coupling constant K and the spread a in the 



oscillator frequencies. Sliiino and Frankowicz 23| by a combination of numerical experi 



ments and analysis obtain an approximate K — a plane phase diagram. References 24 



31 



25| 



32J (i.e., Zi = for all oscillators) 



examine the transition between "amplitude death" 

and collective oscillation, explicitly obtaining analytical results for the boundary in K — a 

space separating death and collective oscillation. 

Matthews et al. [26|,|27|, in addition to presenting an extensive numerical exploration, also 
develop an analytical technique for handling the transition to globally coherent oscillation 
rom phase-incoherent individual oscillation with \zi\ > (as in the Kuramoto transition 



18|- 22|); thus this work was the first to include analysis of the effect of amp 



itude dynamics 
m, ui was the 



on this type of transition. In addition, another important result of Refs. 
numerical discovery that near the boundary in parameter space where the transition to 
collective behavior occurs, this collective behavior can be rather complex, including period 
doubling cascades, chaos, quasiperiodicity and hysteresis. Further, sufficiently far above the 
boundary it was found that steady oscillatory behavior prevails (as in the Kuramoto model). 

Reference 28|] introduces a situation that the authors call "aging" in which there are 
two populations, each described by an equation of the form of (171) (with ai = a — K and 
r = K); the "old" population has dj = —Oo < (corresponding to amplitude death at 
K = 0), and the "young" population has at = ay > 0; Ui was taken to be the same constant 
fl for all old and young oscillators (see also 29| which allows distinct old and young natural 
frequencies, Ui = Qo and Qy); and behavior was investigated as a function of the ratio of the 
populations of old relative to young. In the set up of Refs. [28|, l29[, due to the homogeneity 
of frequencies, the transition problem reduces to the analysis of two coupled Landau-Stuart 
equations. 

Nonlinear behavior of large systems of identical Landau-Stuart oscillators was considered 
by Refs. [37|-|4l|, which highlight the possible occurrence of "clustered states," such that 
oscillators in the same cluster all behaves identically, but Zi{t) ^ Zj{t) if i and j are in 
different clusters. 

The rest of this paper is organized as follows. Section IIIII derives the characteristic 
equation governing linear stability of perturbations from the (z) = state. Section IIVI 
evaluates the characteristic equation for the case of a Lorentzian frequency distribution, 
g{u}) = [7r(l + uj'^)]^^. Section IVl evaluates the effect of spread in a on linear stability in the 
case of Lorentzian g{u}) and no nonlinear frequency shift (7 = 0). Section IVTl evaluates the 



8 



effect of nonlinear frequency shift (7 7^ 0) on stability in the case of Lorentzian g{u}) and no 
spread in a. Section lylTl considers a flat-top frequency distribution, g^u) = f/(l — |w|)/2 
(where f/(») denotes the unit step function) and investigates whether the qualitative behavior 
found in Sees. IVl and IVll is affected by this change in the form of g{uj). Section IVIIII discusses 
behavior above the instability threshold for cases when there is no spread in the nonlinear 
parameters /3 and 7 of ([1]) (as in Eq. ([7])). Section [IX] studies stability of the corresponding 
nonlinear solutions in the limit of large coupling, T/Tc — )■ 00, where Tc denotes the critical 
value of r at which the {z) = state becomes unstable. A primary issue addressed in Sees. 
IVIIII and IIXI is the explanation of why, for sufficiently small 7, macroscopic solutions become 
purely oscillatory with constant amplitude as F/Fc is increased (referred to as the "locked 
state"), and shows that multiple-clustered states with complex dynamics can occur at large 
7. Conclusions and further discussions are given in Sec. IXII 

III. LINEAR STABILITY OF THE {z) = STATE 

We consider Eqs. ([7]) in the limit A^ — )■ 00. In this case, there is a solution corresponding 
to zero value of the order parameter {z). For (z) = 0, Eq. fITaj) has the solutions 

Zi = for ai < 0, (10a) 

Zi = ^/aiexp{i[{ui — 'yai)t + 6'oj]} for ctj > 0. (10b) 



We express the order parameter (z) as 



(z) = {z)^ + {z)+, (11a) 

)- = ^ E ^- (lib) 



i,ai<0 



(^)+ = ^ E ^- (11^) 



N 

i,ai>0 



That is, {z)_ and (z)^ denote the contribution to the order parameter from oscillators with 
ttj < and ai > 0, respectively. Note that (llOap implies {z)_ = 0, while ( llObj) implies 
{z)+ = if A^ — 7- 00 and the angles 6'oj are uniformly distributed in [0, 27r]. Thus, by (I Hal) . 



we see that {z) = is indeed a self-consistent solution of the system ([7]) for A^ — )■ oo. We now 
ask whether this solution is stable to small perturbations. If it is not, then the state {z) = 
will not persist, and global collective behavior will result. We denote the perturbation of 
the order parameter by 



{6z) = {6z). + {6z)+, (12a) 

where Szi is a perturbation from the unperturbed orbit dynamics given by Eqs. (TTUj) . 

Calculation of {6z)^. Considering oscillator i for which a^ < 0, and perturbing Eq. ( 17a|l 
about Zi = 0, we obtain the linearized equation, 

-—^ = {ai + ioJi)6zi + r{6z). (13) 

Assuming exponential time dependence of the orbit perturbations, Szi ~ exp(st), Eq. ( fT3|l 
yields 

Sz, = . '^}^^^ . . for a, < 0. (14) 

[s + \ai\ — too) 

Thus 

{6z)_ = r{6z) / / ^ II . dadu, (15) 

J-ocJ-oo {S+ l"l -«W) 

where for A^ — )■ oo we have replaced the sum over i in ( jllbl) by integration over w and a 
weighted by the pdf's g{co) and h{a) [Eq. ([6])]. Note that the a integration in ( 1T5l) runs 
from a = — ooto« = and thus includes only those oscillators for which a < 0. 

Formulation for calculating {Sz)+. We begin by re-expressing Eq. (17al) in polar form, 
Zi = piexp{i9i) where pi{t) and 9i{t) are real. 



^ = aip, - p3 + rRe{e-^'{z)}, (16a) 

-^ = a;,-7P? + -/m{e-^(^)}, (16b) 
dt Pi 

{z) = {pe'). (16c) 
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We now introduce a pdf for the state variables (p, 6) and parameters {u, a) which we denote 
by 



f{p,9,u,a,t). 



Thus 



27r /"OO 



JO 



fdpdO = g{u)h{a). 



By conservation of the number of oscillators and Eqs. ( TT6l) , / satisfies the following continuity 
equation, 



df 


+ 


d f 


dt 


dp\ 






3 f 




+ 


de \ 



ap-p' + -{e-^'{z) + e^'{zY) 



f 



UJ 



^p' + ±-{e-^\z)-e^\zY) 



(17) 



/ =0, 



where 

{z) = III pe'^fdpdOduda. 

For a > the time independent incoherent ((;z)+ = 0) solution of ( fTTl) and ( fTSj) is 

g{uj)h{a) 



/o 



277 



-5{p- ^/a). 



(18) 



(19) 



We now introduce a perturbation to the solution flTOl) . 



/ = /o + e^*-^''(5/ + {O.P.r.}, 



(20) 



where {O.P.T.} denotes "other perturbation terms" whose 6 variation is proportional to 
exp(m^) with n ^ —1. These other terms do not contribute to (z)+ [see Eq. ( ITSl) ] and 
so turn out to be of no consequence to what follows. Inserting ( 120|) and ( IT9l) into ( IT71) we 
obtain for a > 



S — 2CJ 



r(fe)^(a;)/i(ct ) [5(p-v/a) 
An 



6'{p- 



(21) 
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where 

5'{p - y/a) = —5{p - v/«). 

Calculation of {Sz) + . We now solve (I2T1) for 6f. To do this we assume a solution of the 
form 

Sf = ^^^^^^^^^^[ci(a;, s)5ip - V^) + c,{uj, s)5'{p - v^)], (22) 

and substitute this assumed form into fl2T|) . Using the delta function identities 



F{p)5{p - v^) = F{^)5{p - v/^), 

F(p)5'(p - v^) = F(v^)5'(p - v^) - F'(v^)5(p - v^), 

(where the second of these identities follows from differentiating the first), Eq. fl2T|) yields 

Ti + T2 = ^5 - 5' (23) 

'a 



where 5 = 6{p — i/a); ^' = ^'(p ~ v^)? ^1 results from the first term on the left hand side 

of (EH), 



Ti = {s - iuj +Z7p^)(ci(5 + C2(5') 

= (s — zw +z7Q;)(ci(5 + C26') — 2z7A/ac2(5, 



and To results from the second term on the left hand side of (EI 



T2 = |-{(a - p^)p{ciS + C2S')} = 2ac2S'. 
op 

Separately equating coefficients of 6 and 6' on the two sides of ( 123|) . we obtain two linear 
equations for the coefficients Ci and C2. Solution of these equations yields 



1 s — iuj + 2a — i'ya 

Cl 



^/a s — iuj + 2a + i'ya s — ioo + i'ya ' 
1 
s — iu + i'ya + 2a 

Insertion of (|22|) with these expressions for ci and C2 into (IT8|) then yields {6z)^, 

12 



/ + 00 
dug{u) 

(24) 
r+°° (s ~ iu + a)h{a) 1 ^ ^ 

/o [s — iu + 2a + i'ya] [s — iu + i'ya] J 

Note that the a integration in (^^ is only over positive a (i.e., the integration runs from 
a = to a = cxD.) 

Equation determining s. Inserting flTSj) and fl24l) into (I12ap we obtain, 



_;^ / f^^ [s — iu + a)g{u)h{a)dadu 



ooJo [s-iuj + 2a + i-fa] [s - iu + 27a] 

J-ooJ^oo s + \a\-iu 
By causahty, this expression for the dispersion function D{s), as well as our previous results, 
Eqs. (fT5|) and (!24l) . for {Sz)- and ((^-2)+, are defined with Re{s) > 0. This implies the u- 
integration contour should pass above all poles in the complex cj-plane. We note that for 
Re{s) > the w-integration poles in (I15p . (^^ and (12 5 p all lie in the lower half w-plane. Since 
we are interested in the occurrence of instability, and instability corresponds to Re{s) > 0, 
the form giving D{s) by ( l25l) is sufficient for our purposes {D{s) for Re{s) < can be 
obtained by analytic continuation, from the Re{s) > result). 

IV. LORENTZIAN FREQUENCY DISTRIBUTION 

As discussed in Sec. [H and as we will verify by the example in Sec. I VIII we believe that 
different monotonically decreasing, continuous frequency distribution functions g{uj) often 
(but not always Ref. [33|) yield similar qualitative behavior, and we, therefore, specialize 
here to one such g{uj) that allows easy analytic evaluation of the integrals over u, namely, 
the case of Lorentzian g{u}), 

^(^) = 1^1 = 2^ {^ - ^} ' ('^) 

where we have adopted a normalization of t, T and a so that the half-width of g{u}) is one 
{g{0) = 2g{l)). Since Re{s) > 0, the only w-pole of the integrands in (l25l) that is located 
in Im{uj) > is the one at a; = i [see Eq. (|26|) ]. In addition, the magnitudes of the 
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integrands behave like \uj\~^ for large l^l. Thus, we can deform the w-integration path by 
shifting it upward into the complex w-plane, letting Im{u) along the path approach +00. 
The integration then yields the residue of the pole at u = i, 

/•+°° {s + 1 + a)h{a)da 

~ Jo [s + 1 + 2a + i^a][s + 1 + t^a] 
n[a)aa 



J s + |q;| + 1 

In Sec. |V] we investigate conditions under which Eq. (!27|) predicts instability (i.e., 
existence of a solution to D{s) = T'^ with Re{s) > 0). 

V. CONDITION FOR INSTABILITY OF THE {z) = STATE: THE EFFECT OF 
A SPREAD IN THE GROWTH RATES IN a 

In this section, we consider the case where there is no nonlinear frequency shift (i.e., 
7 = 0), with g{co) being Lorentzian. Using a generalization of the technique in Ref. [2a] (see 
proof of their Theorem 2), it can be shown that the solution of D{s) = l/F is real. Thus, 
as we pass from stability to instability, s goes through s = 0. This results in the following 
general condition for instability, 

r > ^T, (28) 

and (!27ll and (1281) imply that instability occurs when F exceeds the critical value Tc given 
by 



c 



" ^ h{a)da + / h{a) \ . (29) 



/o 2a + 1 J _ao I — a 

As a simple reference case, we first consider (!29|) when there is no dispersion in a 



h{a) = 6{a — a) 



in which case we obtain 



^ ^ (2a + l)/(a+l), fora>0, 
1 + |a|, for a < 0. 

The resulting phase diagram is given by the black line in Fig. [21 This result (with the 
different parametrization used in Eq. ([H])) has been previously obtained in Refs. |26|, l27l |. 
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Note that Fc — ?■ 2 as « — )■ +00. The value Fc = 2 is the critical coupling value that applies 
for the Kuramoto model with a Lorentzian frequency pdf, Eq. fl26l) . The applicability 
of the Kuramoto result for large a can be understood from Eq. fll6al) with F neglected, 
dp/dt = ap — p^, which when linearized about the incoherent equilibrium value, p = a/^, 
yields d6p/dt = —2a6p. Thus perturbations from p = a/^ relax at the exponential rate 2a, 
and, for large a, this rate becomes much faster than the other relevant time scale, namely, 
the spread in u (which we have normalized to 1). Hence, for a ^ 1, the oscillator amplitude 
is essentially frozen, and the Kuramoto oscillator description is valid. As shown in Fig. [2] 
and Eq. (150]) . when a ^ 1 does not hold, the effect of amplitude dynamics is to reduce 
Fc (for a > 0) from the Kuramoto value with the reduction increasing to a factor of 2 as 
a — )■ 0"'' (Fc = 2 at a — )■ +00 in comparison with Fc = 1 at a = 0). An additional interesting 
point is that comparison of the black line in Fig. [2] with the results for the phase diagram 
in the case of Belousov-Zhabotinsky system of Ref. [^ (see discussion in Sect. [T]) shows a 
striking qualitative similarity between the two (e.g., see Fig. 3 of Ref. |2|). 

Referring to Eq. (1301) and Fig. [2], we see that there is a sharp transition in behavior as 
a crosses a = 0. In particular, the (z) = state for a < results from the fact that Zi = 
for all oscillators, while for a > all oscillators have \zi\ = y/^ > and (z) = results 
from incoherence of the individual oscillator phases. This sharp transition in behavior is 
reflected by the discontinuity of the derivative, dTc/da, at a = 0. The sharp nature of the 
transition at a = is, however, a nonphysical artifact of the assumption of no dispersion in 
the individual oscillator growth / damping rates used in obtaining fl5D]) . In typical physical 



I2 



situations, such as the experiment in Ref. J2| (see discussion in Sect. HI]), dispersion in a is 
to be expected (Fig. [1]). To simply illustrate its effect we consider the example where h{a) 
is uniform within some range 6a about an average value a, 

h{a) = {2Sa)-^U{Sa-\a~a\), (31) 

where U{x) denotes the unit step function; U{x) = 1 for a; > and U{x) = for x < 0. 
Using f ET]) in f l2^ . we get for a > 6a, 



F-i 



2Sa 
for a < —6a, 



1, fa + 6a + l/2 
da -\ — In ' 



4 \a- 6a +1/2 



(32) 



^-2LKl^f^l^ (-) 
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a 



FIG. 2: (Color online) The effect of dispersion in a (7 = 0) with a Lorentzian g{uj). Stabil- 
ity/Instability regions of F — a space for several different values of spread 5a in the linear growth 
parameter a with mean a [Legend: Black line {6a = 0), red line {6a = 0.5), blue line {6a = 2)]. 



and for |a| < 5a, 



25a 



a -\- 5a 1 

h ln(l -a + 5a) + - ln(l + 2a + 25a) 



(34) 



As the dispersion in a, 5a, approaches zero, (1311) becomes a delta function, and Eqs. (132|1 - 
( !34l) reduce to ( !30l) . The other two Hues in Fig. [2] show the phase diagram from Eqs. 
(!32l) -( l34l) for two more values of 5a. For 5a > Q the discontinuity in dTc/da (which occurs 
for 5a = a.t a = 0) is removed by dispersion in a, and the sharp transition that occurs at 
a = (black line in Fig. [2]) is now smoothed out 3J]. Further, it is also noticed that the 
minimum of F^ rises and shifts from a = when 5a = to a > when 5a > 0. 



VI. CONDITION FOR INSTABILITY OF THE (z) = STATE: THE EFFECT OF 
A NONLINEAR FREQUENCY SHIFT 

We now address the effect of nonlinear frequency shift, 7 7^ 0, and we consider the simple 
case of no dispersion in a, h{a) = 5{a — a) again for the case of Lorentzian g{uj). We note 
from Eq. ([7]), if the distribution of Ui values is symmetric, then positive and negative values 
of 7 are equivalent {zi — )■ z*). Thus, we consider 7 > only. As is evident from Eq. fl27|) . 
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7 has no effect on the hnear theory for a < 0, and, consequently, the result for Tc given by 
Eq. fl30|) still applies for a < 0. For a > 0, however, the effect of a nonlinear frequency shift 
can be substantial. Equation ( HT} for h{a) = 6{a — a), a > gives 

r-^ = Dis) = , ^±1±^ ,, (35) 

which yields a quadratic equation for s, solution of which can be used to obtain stability 
boundary curves in F — a space. At the transition point, Re{s) goes through zero. Substi- 
tuting s = iQ into Eq. fl35|) and separating the real and imaginary parts, Tc and Q are then 
given by the solution of the following pair of equations 

1 - fi2 + 2a(l - 7^) - aV = ^c{l + «), (36a) 

2{l + a){n + a^)=T^n. (36b) 



When 7 = (note fi = for this case), the solution for the critical coupling strength of fl36|) 
is given by 



r. = ii^, (37) 

1 + a 
by which f l36ap can be rearranged to give 

Fc = Fco -— — , (38) 

1 + a 

which shows that the effect of 7 is always to decrease Fc. Figure [3] shows the values of F^ 

as a function of a for several different values of 7 (7 = plotted in black, 7 = 2 plotted in 

red, and 7 = 4 plotted in blue). By solving for Q in (]36bp and substituting it back in (I36ap . 

we obtain 

a'yTc 



cO 



1 + a 



Fc-2(l + a) 



(39) 



Equation (1391) shows that Fc — )■ 2 as a; — t- +00. As seen in Fig. |3l increasing 7 eventually 
moves the minimum of Fc below one and shifts the location of the minimum into a > 0. 

VII. THE EFFECT OF THE FREQUENCY DISTRIBUTION FUNCTION 

In Sees. |V] and |Vl] we consider the effect of a spread in a and of a nonlinear frequency 
shift for the illustrative case of a Lorentzian distribution function of the oscillator natural 
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FIG. 3: (Color online) Stability / Instability regions of F — q space for several different values of 
the nonlinear frequency shift parameter 7 [Legend: Black line (7 = 0), red line (7 = 2), blue line 
(7 = 4)]. Notice that the three lines coincide when q < 0. 

frequencies, Eq. ( 126|1 . We now ask how might these results be altered if a different frequency 
distribution were used. We note that the Lorentzian decays rather slowly for large u, g{uj) ~ 
u^^. Thus to test dependence on the form of g{uj), we will examine another distribution 
function which is very different from the Lorentzian, in that it has a sharp cutoff to g{uj) = 
as u increases. In particular, we will consider a "fiat-top" distribution, that is uniform in 
— 1 < w < 1 and zero otherwise, 



g{u) = -Uil-\u\), 



(40) 



where U{x) is the unit step function. In spite the qualitatively different large |a;| behavior 
of the Lorentzian and the fiat-top g{uj) distributions, we will find that the resulting stability 
conditions show qualitatively similar behavior. 

The calculation of Tc with g{ijj) given by fj40|) is done by using fl25|) (see the Appendix). 
In Fig. IHwe show the dependence of Fc on a for several different values of Sa, where h{a) 
is given by ( 13T]) and 7 = for all oscillators. A comparison between Fig. H] and Fig. [2] 
reveals remarkably similar dependence, apart from a difference in the vertical scale due to 
different functional dependence of g{u) [35|. Next, we consider the dependence of Fc on 7 
when h{a) = 6{a — a) with g{u) given by (j40|) . Figure [5] shows the dependence of Fc on a 
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FIG. 4: (Color online) The effect of dispersion in a with a uniformly distributed g{u;), Eq. (|40p . 
Stability/Instability regions of T — a space for several different values of spread 5a in the linear 
growth parameter a with mean a [Legend: Black line {6a = 0), red line {6a = 0.5), blue line 
{6a = 2.0)]. 

for several different values of 7. The black line shows the result when 7 = 0, which is the 
same black line in Fig. HI The other two lines are obtained by numerically solving Eq. (1A7I) 
in the Appendix when 7 7^ 0. In comparison with Fig. [3], we see similar dependence in that, 
as 7 increases, Tc decreases. 

VIII. NONLINEAR PHENOMENA ABOVE THE INSTABILITY THRESHOLD 
WITH FINITE a-SPREAD AND NONLINEAR FREQUENCY SHIFT 

In the previous sections, we calculated the critical coupling strength F^ marking the 
onset of instability of the quiescent state {{z) = 0). Above the critical value F^, we find 
that Landau-Stuart oscillator networks exhibit a rich variety of collective behavior. We now 
briefly review past work on the nonlinear behavior found above the instability threshold. 

Matthews et al. [27i] studied nonlinear collective behavior in the special case that aj = 1 
and 7j = for all oscillators {j = 1, 2, ■ ■ ■ A^), and g{uj) takes on several different functional 
forms. An important observation in that paper is that the system behavior can be quite 
complicated for a range of F not too far above T^. For example, they found period doubling 
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FIG. 5: (Color online) Stability / Instability regions of F — q space for several different values of 
the nonlinear frequency shift parameter 7 [Legend: Black line (7 = 0), red line (7 = 2), blue line 
(7 = 4)]. Notice that the three lines coincide when q < 0. 

cascades to chaos, large amplitude oscillations, quasiperiodicity, and hysteretic behavior 
close to the boundaries between different macroscopic states. However, when F is sufficiently 
far from Fc, the system was always found to settle into a steady oscillatory state, {z) = 
constant x exp(if2t) for some constant Q. We refer to this as a "locked state," which we 
define as a solution of ( ITTll and ( ITSl) for which the oscillator distribution / has dependence 
on (p, 9, t) of the form / = /(p, 6 — fit); i.e., the entire distribution rigidly rotates about the 
origin of the complex 2;-plane with the uniform rotation rate fl. 
When the nonlinear frequency shift parameter 7 is nonzero 



37 



401 1 , the system can exhibit 



additional types of complicated coherent behavior. For example, Refs. |37H4l| studied 
systems closely related to Eq. (JTj), but with homogeneous parameters. An important feature 
found in those references is the tendency for the system to form clusters (a "cluster" in this 
case is defined as a group of oscillators which behave identically). Further, depending on 
parameter values and on initial conditions, the systems can form cluster states of varying 



sizes. In Refs. 



and 



the authors also found chaotic behavior. 



We emphasize the finding of Ref. 271] that, for zero nonlinear frequency shift 7j = the 
system always goes to a locked state attractor when F is sufficiently large. Consistent with 
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this, we find that when we include spreads in a, and w, and simultaneously allow f3j = ^ j^ 
and 7j = 7 7^ 0, it is the case that, as F is increased, there is always a locked state that the 
system may settle into. Furthermore, we analytically prove that this locked state is the only 
large F attractor (as in [27] which has 7j = 0) provided that the nonlinear frequency shift 7 
is not too large, but we also find that other coexisting attractors may be present if 7 is large 
enough. This will be discussed further in Sees. |IX] and |Xl As an example of a locked state. 
Fig. [6] shows snapshots of the long-time asymptotic distributions of the oscillator z— values 
obtained from numerical simulations of Eq. ([7]) with g{uj) given by (HOl) . A^ = 50000, h{a) 
given by (jH]), a = 0.5, 6a = 1.0, 7 = 0.5 (corresponding to Fc = 0.89), for successively 



larger values of F/F^, all of which are large enough that a locked state is achieved (Fig. 6(a 



F/F, = 2, Fig. [6(b)| F/F^ = 10, Fig. [6(c)| F/F, = 100). Note that, as appropriate for a 
locked state, as time increases, these snapshots rotate uniformly about the origin at a fixed 
angular rate Q. 



We see in Fig. 6(a) that the distribution has finite spreads in both the magnitude and 
phase of z. Examination of the solution shows that oscillators with smaller (larger) natural 
frequencies u tend to occur on the clockwise (counterclockwise) side of the distribution, while 
larger (smaller) a tend to occur at larger (smaller) \z\ for fixed argument of z. Previous 
works (e.g., [23]) did not consider a distribution of a and consequently did not find a spread 
in \z\ at constant argument of z (i.e., the oscillators are distributed along a curve in the 



complex z-plane). Comparing Figs. 6(a), 6(b) and 6(c), we see that the spread in z/\{z)\ 



becomes smaller and smaller as F/Fc increases. In fact, we argue in Sec. IIXI below that 
one of the stationary states of this system is when this spread goes to zero as F/Fc — ?■ 00 



(note the greatly magnified scale for Fig. 6(c)). Note that the oscillators in Fig. [6] are 



contained within a single region, and we subsequently refer to such states as single-cluster 
locked states. 

Figure [3 illustrates an example of the occurrence and evolution of a non- locked dynamical 
attractor at lower F/Fc with other parameters the same as those in Fig. [6l In particular, 
Fig. [7] shows 1(2)1 (top panel) and Re{z) (bottom panel) versus time, after the system has 
settled into an attractor for F/Fc = 1-07. (Note that for a locked state, \{z)\ is constant, 
and Re(2;) varies sinusoidally in time.) 
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FIG. 6: (Color online) Locations of 50000 oscillators (black) in locked states with different 
r/Fc- Twenty oscillators (red cross) of parameter values evenly spaced simultaneously in {a,uj) G 
[—0.5,1.5] X [—1,1] are highlighted, i.e., the oscillator with {a,uj) = (-0.5,-1) is located at the 
minimum radius position, and the oscillator with {a,uj) = (1.5,1) is located at the maximum 
radius position, and other oscillators of intermediate parameter values are distributed in between 
{N = 50000, a = 0.5, 6a = 1.0, 7 = 0.5; random initial conditions). 

IX. LARGE COUPLED LANDAU-STUART OSCILLATOR NETWORKS IN THE 
STRONG COUPLING LIMIT: SINGLE-CLUSTER LOCKED STATES 

In what follows, as in all other previous references (except for the weak "coupling limit" 
treatment in Ref. (36|]), we consider the case where there is no spread in the nonlinear 
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FIG. 7: Time evolution of |(2;)| (top panel) and Re(z) (bottom panel) for a system of 500,000 
oscillators (Parameters: a = 0.5, Sa = 1.0, 7 = 0.5, T/Tc = 1.07; random initial conditions.) 

coefficients (7^ = 7 and Pj = 1 for all j). In this section we ask why the single-cluster locked 
state is an attractor for large enough F. 

In order to analytically show that a single-cluster locked state attractor must exist for 
homogeneous nonlinearity parameters 7^- = 7 and f3j = 1 at sufficiently large r/F^, we now 
consider very large F/Fc approximated by taking the limit F/Fc — t- 00. In particular, using 
this limit, we will show the existence of a simple single-cluster locked state and we will 
demonstrate that it is stable. In Sec. |X]we will show that the single-cluster locked state is 
the only attractor of the system if 7 is not too large, but that, when 7 is larger, there can 
be other coexisting attractors of various types composed of multiple clusters. 

When F ^ Cij.Uj for all j, system ([7]) reduces to 



dt 



il + t^)\z,\'z,+T{z). 



(41) 



Here we have assumed that \zj\ ^ 1 in the F — )■ 00 limit. This will be subsequently verified. 
Alternatively, if the Uj are uniform, Uj = u, and F ^ a^, even if F ^ a; does not apply, 
we can still obtain Eq. fj4T]) via elimination of u through the transformation Zj — )■ ZjC"^^^. 
Thus, in this limit, the dynamics is determined by the coupling to other oscillators and the 
nonlinear characteristics of the individual oscillators, rather than by the linear properties 
of the individual oscillators. This is consistent with our numerical tests in Fig. |6l which 
suggests that as F/Fc — ^ 00, the distribution of oscillators approaches that of a system of 
homogeneous parameter values, with the effects of spreads due to a and u going away (see 
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Fig. 6(c)). We now divide Eq.( BT|) by F, and redefine variables as 



% = ^, (42a) 

i = Ft. (42b) 



Thus, each term in Eq. ( BTj) scales as F'^/^ justifying the neglect of the other terms in Eq. 
d?]). Equations (HTj) become 

N 



dz- 1 



fc=l 

Making the single-cluster locked state ansatz Zj = p{t)exp[i9{t)] gives 

^ = pil-p% (44a) 

§ = -IP'- (44b) 

at 

This yields the time asymptotic attracting solution, 

io(t) = e-'~^K (45) 

To analyze the stability of (H5ll . we perturb ;zo to zq + e~^^^5zj. From ( H3ll . the dynamics 
of (5zj is governed by 

-^(5% = (-1 - ii) {5z, + ^*) - 5z, + j^Y. ^^^^- (46) 

A; 

where * denotes complex conjugation. Similarly, we have 

jM] = (-1 + ^7) {^h + <55;) - bz] + 1 ^ 5i^ (47) 

Equations (146!) and (H71) can be regarded as two independent equations for bzj and bz* 
respectively. To study the stability properties of bz^ and (5z*, consider bzj ~ bZj{s)d^^ and 
bz] ~ (5i'*(s)e'*, for which Eqs. (36]) and (iTD give 

(s + 2 + i7) (5Zj- + (1 + 27) bZ* - {bZ) = 0, (48a) 

(s + 2 - ^7) 5Z* + (1 - ^7) 5Zj- - (5Z*) = 0, (48b) 
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where {bZ) = N ^ J2k ^^k, {SZ*) = N ^ J2k ^^l- Summing over j we obtain 



T 



(SZ) 
{6Z*) 



(49) 



where 



T 



(50) 



(s + 1 + Z7) (1 + i^) 
(1 - 27) (s + 1 - i^) 

Equation (^ imphes that either (i) detT = 0, or iu){6Z) = {5Z*) = 0. Possibihty (i) 
gives s{s + 2) = 0, yielding s = and s = —2. Physically, the neutrally stable root, s = 0, 
corresponds to a uniform, rigid rotation of the phases of all the Zj. If possibility (ii) applies, 
Eqs. (I48ap and fl48bp become 



^ 



SZi 



6z; 



0, 



(51) 



^ 



(52) 



where 

(s + 2 + n) (1 + n) 

(1 - i7) (s + 2 - i^) 

Since det\E' = (s + l)(s + 3), we obtain two additional roots s = —1 and s = —3. Because 
the allowed perturbations in case (ii) are restricted to lie in the 2(A^ — 1) dimensional space 
specified by the two constraints, {6Z) = and {6Z*) = 0, the multiplicity of each of the 
roots s = — 1 and s = —3 is A^ — 1. Since there is no root with Re{s) > 0, the equilibrium 
is stable. Hence, the single-cluster locked state is stable. 



X. LARGE COUPLED LANDAU-STUART OSCILLATOR NETWORKS IN THE 
STRONG COUPLING LIMIT: CLUSTER STATES 



A. Regime of global attraction for the single-cluster locked state 



In Sec. IVIIII we numerically suggest the tendency of system (H3ll to form a single-cluster 
locked state when F is sufficiently large, and in Sec. llXI we have shown that such an attractor 
always exists at large F. In this section, we give a sufficient condition for this state to be 
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the only attractor of the system. Consider any two oscillators, m and n, in ( H3|) . Let 

2 ' 



(53a) 
(53b) 



Then the dynamical equation for the separation between the two oscillators, 6, can be 
immediately derived from Eq. ( l43ll as 

A^ = -(1 + Z7) {2\I\^S + 6*F + \6\^6) . (54) 

Letting 5 = \5\e^'^'^ and z = \z\e^'"^ , substitution into (l5l|) yields 

^\5\ = -|l|Vl[2 + cos(z.3) +7sin(z.3)] - |5|^ (55) 

at 

where v^ = 2(z/i — 1/2). The trigonometric terms on the right hand side of Eq. (!55|l can be 
combined, giving 

d 



-.\6\ = -|l|Vl[2 + Vl + 7'cos(i/4 - ^^3)] - |<5|', (56) 

at 



where z/4 is determined by the conditions cos(z/4) = l/yl + 7^ and sin(z/4) = 7/vl + 7^. 



From Eq. fl56|) . we see that (i|5|/(it < if 2 > a/1 + 7^, or equivalently 7 < y^. Thus, 
if 7 < v3, the attractor of system (l43ll must occur as a single-cluster, and, as shown in 
Sec. \IX.\ a single-cluster attractor must be a locked state ( l45l) . However, we shall soon see 
that if 7 > v3, then Eq. ( H3l) has the possibility of attracting solutions other than the 
single-cluster locked state. This technique [Eqs. ( 153|) - (!56|) ] can also be employed for related 
problems such as when the coupling is complex, F — )■ Fe*^, or when the coupling is in the 
form in Eq. ([8]). 

B. Cluster States 

We now wish to investigate the possible existence of attracting states for (1431) composed 
of C > 2 clusters, where each cluster is labeled by a subscript c = 1, 2, ■ ■ ■ ,C. Each cluster 
c has Nc oscillators in identical states Zc, such that if oscillators i and j are in cluster c, then 
f j = Zj = Zc, and Ni + N2 + ■ ■ ■ Nc = N. Letting C,c = N^/N be the fraction of oscillators 
in cluster c, we have that 
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{Z)=irZr + --- + icZc. (57) 

and that Eq. ( l43l) reduces to C equations for the C complex cluster variables Z^ (c = 
1,2,---,C), 

dZc 



dt 



l + t^)\ZfZ, + J2^cZc. (58) 



c=l 

Two questions pertaining to such cluster states are (i) what are the attractors of (ISSj) . and 
(ii) given an attractor of flSSl) . are the clusters internally stable? Question (ii) asks whether, 
if we consider oscillators in cluster c and individually independently perturb each of them 
from their common value Zc, do they relax back to a common value? (This question was 
considered for the one-cluster locked state in Sec. HXl ) 

The question of the existence of cluster state attractors is a general one applicable to any 
large system of identical dynamical units that are coupled via a global field (e.g., in our case 
(Z)). In particular, this type of consideration was introduced by Kaneko who considered 



coupled maps |42l. |43|. 



In our numerical experiments we have always found that, at long time, the solutions of 
(l43ll settles into a finite number of clusters Zc{t). We caution that this does not necessarily 
mean that attractors of fH5]) always occur in clusters, but rather that we have so far not found 



non-clustered long-time states. References [37|-|40| consider a globally coupled Landau- 



Stuart system with homogeneous parameters across all oscillators (our (143!) is a special case) 
and find both clustered state attractors and "scattered state" attractors, where by "scattered 
states" they mean that, at any given time, no two oscillators have exactly the same state. 
We, however, have not seen scattered state attractors, and we conjecture that they do not 
exist for (H5]) . Along these lines, we now show a partial result implying that scattered states 
cannot have scattering that is over an area in the 2;-plane. That is, in the limit N ^ oo the 
distribution function / appearing in flT7|) must be singular in the sense that it is concentrated 
on a set of zero Lebesgue measure in 2;-space, equivalently (p, 6) space. Examples of zero 
Lebesgue measure sets are a set of distinct points (like our clusters), a curve, or a fractal 
set of dimension between one and two. Indeed the scattered states seen in the figures of the 
p..evious .efe..e,.es H-fl (e.,. Fig. 1 of Ref. Q) appea.. to ou. eye to be eithe. .acta. 



distributions with dimension near one or distributed along a convoluted curve (based on 

27 



Ref. J4J], we suspect that the first of these alternatives apphes). The demonstration that / 
must be singular follows simply from f lT7|) by introducing F = f/p and rewriting flT7|) in the 
form 

where p and 6 are the two quantities in flT7|) appearing within the square brackets with 
the linear oscillator parameters, a and u, set to zero [to correspond to fH3|) ]. Note too 
that F = F{p,9,t); in particular, unlike our more general set-up, Eq. (TT7|) . F does not 
incorporate distribution in parameters as we have fixed a and u at zero. According to 
(!59|) . following the characteristics, dp/dt = p, dO/dt = 6, of the partial differential equation 
f lT7|) with a = cj = 0, F increases continually at the exponential rate 4p^. Thus, since 
J fdpdO = J FpdpdO = 1, we immediately conclude that in the long time limit, F and hence 
/ must concentrate on a set of zero area (zero Lebesgue measure) in (p, 9) space. Proof 
that our system fH3|) does or does not always yield cluster state attractors remains an open 
problem. 

C. Two-cluster states 



For C = 2, Eq. ([58]) yields 

dZi _ 

-(l + z7)|Z2pZ2 + 6^i + 6^2. (60b) 



„ -(l + z7)|Zi|^Zi + ei^i + 6^2, (60a) 

dt 

dZ2 

di 



Letting Zi = pie^^^ and Z2 = P2e*^^, and defining the relative phase difference = 6*1 — 62, 
Eq. ( 160|) yields three real equations, 



-JT = Pd^i- Pi)+^2p2Cos{(j)), (61a) 

dt 

-j^ = p2{^2-pl)+^ipicos{(j)), (61b) 

dt 



-jj = -lipl - Pi) - sin(0) 
dt 
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. Pi P2. 



(61c) 



Note that pi = p2 = 1 and = is a solution of these equations and corresponds to the 
single-cluster locked state. We want to obtain two-cluster solutions {Zi ^ Z-i). Although 
we cannot rule out chaotic or two- frequency quasiperiodic two-cluster solutions of Eq. ( 16T|) . 
so far our numerical investigations of Eqs. ( 1611) have only found fixed point attractors 
(i.e., uniformly rotating two-cluster locked states) and periodic attractors. In the next two 
subsections we discuss the fixed point solutions (Sec. IXC 1|) and the periodic solutions (Sec. 
IXC21) . 



1. Two-cluster locked states (fixed point solutions of Eq. /161\} ) 



When the time derivatives in Eq. fl6T|) are all zero, the solutions give locked two-cluster 
fixed point solutions. By setting d/dt = 0, x = pf and y = p\, elimination of cos((/)) between 
Eqs. (16 lap and ( I61bl) gives 

ei(x-|)'-6(y-|)=^(e?-e.^). (62) 

On the other hand, addition of Eqs. (16 lap and (]61bp . and subsequent elimination of the 
trigonometric factor with that in Eq. (I61cp by the identity sin^(0) + cos^(0) = 1 gives 

{^ + y-^f + i\y~^?={i2^- + ii^^ . (63) 

Two-cluster fixed point solutions are given by the intersecting points of Eqs. ( 162|) and ( l63l) . 
An example shown in Fig. [8] corresponding to the parameters ^i = 0.9 and 7 = 4.2. There 
are altogether four intersecting points, but two of them, namely (0,0) and (1,1), do not 
correspond to the answers we seek ((0,0) is the unstable incoherent state and (1, 1) is the 
one-cluster locked state solution). For the other two solutions, indicated as A and B in 
Fig. [HI we find that B corresponds to an unstable fixed point, while A is stable, and our 
numerical tests on Eqs. f l6T]) and f H3|) show that A satisfies both types of stability, (i) and 
(ii) mentioned at the end of Sec. IX B[ Thus A is an attractor. 
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(a) 



0.0006 



(b) 



FIG. 8: (Color online) (a) Two-cluster locked state solutions. The blue line shows the curve for 
Eq. (|63p. and the red line shows the curve for Eq. (j62|) : parameters: ^i = 0.9, 7 = 4.2. (b) Blowup 
of the figure in (a) around (0, 0). 



2. Two-cluster periodic solutions of Eq. ^61\ 



In general, the two-cluster periodic solutions of Eq. f ET]) are hard to obtain analytically. 
However, in the case of large 7, we can proceed using perturbation theory. If 7 ^ 1 and 
pI — P2 is not small, then d(j)/dt is very large, so cos(0) is rapidly varying. Thus, to lowest 
order in 7"^, we can neglect the cos(0) terms in the equations for dpi/dt and dp2/dt, 



dpi,2 
di 



~M'\~S^) 



(6.2 - pZL kl 



(64) 



;^(o) 



and pY' and P2 , are attracted to p^ = a/^, p2 = y/^, respectively. Thus we have 



-^ ~ -7(6 - 6) 7^^ sm(0), 

dt \/4i42 



which, for large 7 has the solution, 



^0 - 7(6 - 6)^' 



vC6 



cos[0o -7(6 -6)^1 • 



(65) 



(66) 



;;(0) , ;;W 



;;(!) 



To next order in 7 , we write pi^2 = Pi 2 + Pi 2 "where the perturbation p\ to the lowest 
order quantity p[ satisfies 



di 



-SeiPr+d^ cos[0o-7to-6)t^■ 



(67) 
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FIG. 9: Simulation study of a two-cluster periodic state; parameters: N = 50000, ^i = 0.8, 7 = 7.2. 

The homogeneous solution of the above equation decays as e~^^^*, and thus does not con- 
tribute to the time periodic attractor. For large 7 the inhomogeneous term varies with 
a period 27r/[7(^i — ^2)] which, since 7 is large, is much shorter than the damping time 
(2(^i)~^. Thus for calculating the inhomogeneous solution, we may neglect the term 2^ipi . 
This yields 

,3/2 
Piit) = - wT^^^T sin[0o - 7(6 - 6)t], (68) 



7(6-6) 



and similarly 



^3/2 

P2^(t) = - -(, ,^ si^[<^o - ^(^1 - ^2)^1- (69) 

7(6 - 6) 

Thus p\2 are indeed small compared to p\2, if 6 ¥" 6 and 7 is sufficiently large. Hence 
we obtain a two cluster time periodic state whose frequency is, to lowest order, 7(^1 — ^2)- 
Figure shows long-time results of a simulation of Eq. fH5]) with 50000 oscillators, and with 
parameters ^1 = 0.8 and 7 = 7.2, for a periodic attractor. Comparison of these results with 
the approximate analytical solution above shows good agreement. 

In order to see why this solution represents an attractor of the full A^ dimensional system 
(H3|) . we consider its Lyapunov exponents. To lowest order the individual clusters are un- 
coupled locked states. We have already shown (Sec. HXl) that a single-cluster locked state of 
a system of A^ oscillators has A^ — 1 negative exponents (having possible values —1, —2, —3) 
and one zero exponent that corresponds to a rigid rotation in the complex plane of the entire 
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system of A^ oscillators (i.e., Zj — )■ Zj exp(2$) for constant $). Thus, to lowest order in our 
7~^ expansion, there are (A^^i — 1) + (A^2 — 1) = A^ — 2 negative Lyapunov exponents and 
two zero Lyapunov exponents. Now introducing small coupling between the clusters (i.e., 
finite 7), the negative exponents will be slightly perturbed by an amount 7^^ ^ 1 and hence 
will remain negative. The only danger of instability is that one of the two zero exponents 
might be perturbed to be a positive number of order 7"^. However, this cannot be the case, 
because the full system must have two zero exponents, and thus the two zero exponents of 
the lowest order uncoupled approximation are preserved. To see this, we note that there is 
one zero exponent corresponding to a rigid rotation of the entire system of A^ = A^^i + A^2 
oscillators. Note that this zero exponent is not present in the three ODE's, Eq. (I6ip . for 
Pi, P2 and (p, since a rigid rotation (6'i — )■ 6'i + $, ^2 — ^ ^2 + ^) does not change cj) = 9i — 62- 
Another zero exponent results from the fact that the time periodic flow, Eq. fl6T]) . has a zero 
exponent corresponding to displacement along its orbit. Thus we conclude that our large-7, 
two-cluster, states are attractors. 

D. Cluster-states with C > 3 

1. Do locked state attractors with three or more clusters exist? 

The above implies that, at large F, we can have both two-cluster and single-cluster locked 
state attractors. A natural question is whether large-F locked state attractors with C > 2 
clusters are possible. We now give a partial answer to this question by ruling out the 
possibility of locked states composed of more than three clusters. To rule out C > 3 locked 
state solutions of Eq. fH3|) . we substitute the locked-state ansatz Z^it) = p^e^^'^'^ e~^^^ into 
( H3|) . where Pc> and 9co are time independent real constants. This yields 

{Z{i)) = {Z{0))e-'^\ and (70a) 

[{l + z^)pl~zn]p,e^'-^ = {Z{0)) (70b) 

Multiplying Eq. ( ]70bl) by its complex conjugate, we obtain 

~pTpt + {rpl-m = \{zm\'. (71) 
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A particular state corresponds to particular values of Q and |(Z'(0))|. Thus Eq. ( !7T|) must 
be satisfied for each individual cluster c composing the state for the same values of fl and 
|(Z(0))|. Since (I7T]) is a cubic equation for pi, there can be at most three real values of 
Pc > 0. Furthermore, Eq. f lTU]) uniquely determines the value of 6co for each value of pc- 
We, therefore, conclude that large-F, locked, cluster states with C > 3 cannot occur. This 
leaves open the question of whether or not three cluster locked state attractors exist. In this 
regard, we note that in our, admittedly limited, series of numerical experiments we have so 
far not seen such attractors. 

2. Periodic, quasiperiodic and chaotic attractors for C > 3 

Considering the C complex ODE's for the C cluster states, Eq. f l58|) . and again intro- 
ducing a polar representation, Zc = Pce^p{iOc), we obtain a 2C — 1 dimensional dynamical 
system of C real equations for dpc/dt and C — 1 real equations for dipc/dt where (pc = 0c — di 
and c = 2, ■ • • C. Again taking 7 ^ 1, we find C — 1 lowest order angle evolutions, 

^' :Aa;„ l^ujc = -i{ic-ii). (72) 



dt 
Assuming that the set of C — 1 frequencies {Ac^c} are irrationally related in the sense that 

the equation, 

c-i 

J]meAwe = 0, (73) 

c=l 

has no solution where the nic are positive or negative integers except for the trivial solution 
where ttIc = for all c, then we can think of the lowest order solution as being (C — 1)- 
quasiperiodic in the 2C—1 cluster-state phase space {pi, ■ ■ ■ , pc] 02, ■ ' ' > 0c}- Application of 
perturbation theory in the small parameter 7"^ is mathematically equivalent to the problem 
of investigating the introduction of small coupling between C — 1 oscillators. 

For example, for three clusters, we have the possibility of two-frequency quasiperiodic 
motion, and the possibly analogous problem of introducing small generic coupling between 



two periodic oscillators was originally addressed by Arnold 45| in his study of the circle map, 
S„+i = {2'KRQ + En + KsinEn) mod 27r, where Rq denotes the rotation number for k = (i?o 
is analogous to A(jJi/Aw2 in our case, and n is analogous to 7"^). Arnold's work resolved the 
problem of the convergence of perturbation theory of coupled oscillators which is plagued by 
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the proliferation of small denominators in higher and higher order terms in the perturbation 



series. Results, both analytical (as by Arnold 45|) and numerical, show that for k < 1, 
attracting quasiperiodic motion continues to exist on a positive Lebesgue measure of the 
parameter space (in our case, the parameter-space is {(.c]!}), but is structurally unstable: 
for any parameter set yielding two-frequency quasiperiodicity, one can find an arbitrarily 
close-by set of parameters where the motion is periodic. Alternatively, one can say that 
two-frequency quasiperiodic attractors exist on a positive Lebesgue measure Cantor set in 
parameter-space, while periodic attractors exist on the complement of this set which is an 
open set (e.g., see |2l| for further discussion). 

Figure dU] shows evidence supporting this scenario. The figure shows the result of numer- 
ical computations of fISS]) for the rotation number defined by 

n. . ,i. MI)::^. (.4, 

versus the rotation number at infinite 7, 

^0 = 1^, (75) 

for C = 3,7 = 30,^3 = 0.2,^1 varied from 0.497 to 0.533, and ^2 = 1-6-6- A 
classic "devil's staircase" pattern is clearly observed, with periodic orbits corresponding to 
frequency-locking plateaus at rational rotation numbers of Re = 2/3, 5/7, 3/4, and 4/5. The 



blowups, shown in Figs. 10(b) and 10(c), reveal further plateaus at Re = 7/10 and 7/9, 
suggesting that (as in Arnold's circle map) plateaus exist for all rational numbers p/q {p, q 
incommensurate integers) with the plateau widths decreasing as q increases. 

The case C — 1 > 3 adds a new ingredient. Again, as 7"^ increases from zero (7 becomes 
finite), (C — l)-quasiperiodic attractors generally persist on a positive Lebesgue measure 
Cantor set in parameter space, but the open complement of this Cantor set now typically 
contains a variety of other types of attractors, including periodic attractors (as for C = 
3}, chaotic attractors, and M-frequency quasiperiodic attractors with M < C — 1 (j46|- 



491]). In particular, the generic existence of structurally stable chaotic attractors resulting 



from perturbations to P-frequency quasiperiodic attractors for P > 3 has been proven by 



Newhouse, Ruelle and Takens 
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Note that all the motions, including the chaotic motions, referred to above occur in 
the context of Eq. (|58|) and are therefore motions of the C clusters. Also note that in our 
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FIG. 10: (a) Rotation number Re from numerical solutions versus the infinite 7 rotation number 
Ro for C = 3,7 = 30,6 = 0.2, and ^i varied from 0.497 to 0.533 with ^2 = 1 - 6 - Cs- (b), (c) 
Blowup of two regions of the figure in (a). 



discussion above of solutions of Eq. ( 135]1 . we have treated the parameters {^c} as continuous. 
While this is formally allowed for Eq. (158|) . we emphasize that if we consider that (l58ll is 
derived from system fH3l) with finite N, then ^ = N^/N can only take on discrete values 
(although they may be very dense for large A^). 
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XI. DISCUSSION AND CONCLUSIONS 

In this paper we have studied some properties of large all-to-all coupled Landau-Stuart 
oscillator networks. The motivation for studying this class of systems is to reveal possible 
generic behaviors of large systems of coupled oscillators where the oscillators have both 
amplitude and phase degrees of freedom. 

In the first half of this paper (Sees. IIimVII|) . motivated by experiments reported in Ref. 



[2|, we studied stability of the incoherent state {{z) = 0) and determined the stability / 
instability boundary for different cases. First, we studied networks with spread in the distri- 
butions of the natural oscillator frequencies u and the linear amplitude growth parameter a, 
but with no nonlinear frequency shift contribution, i.e., 7j = for all oscillators i. Second, 
we studied networks with no spread in the distribution a, but with a constant nonlinear fre- 
quency shift parameter 7 for all oscillators. After establishing a mathematical framework to 
determine the stability / instability boundary, we characterized the changes in the stability 
/ instability boundary that these modifications cause. First, we found that a spread 6a in 
the distribution of a smooths out the discontinuity at a = in the slope of the stability 
/ instability boundary. Second, spread in a causes the minimum of Fc to shift away from 
a = to a > when 6a > 0. Third, increase of the nonlinear frequency shift parameter 7 
monotonically lowers F^. 

Similar to large networks of phase oscillators of the Kuramoto type, large networks of 
Landau-Stuart oscillators with small nonlinear frequency shifts have a tendency to always 
synchronize into a locked state exhibiting steady, constant-amplitude sinusoidal motion when 
the coupling strength is large enough. In order to better understand this behavior, in Sees. 
IVIIIIIXl we studied the limit F/Fc — t- 00. We found that as F/Fc — ?■ 00, Eq. ([7]) reduces to 
Eq. (H3l) , which depends only on coupling among oscillators individually dominated by their 
constant nonlinear characteristics. Considering cluster state attractors of fH3l) we obtained 
the following results. 

1. For sufficiently low values of the nonlinear frequency shift parameter (7 < y/S), there 
is a unique, global attractor that is a single-cluster, locked state attractor. 

2. For larger 7, multiple-cluster attractors can occur, but the single-cluster, locked state 
attractor continues to exist. 



36 



3. For larger 7, two-cluster locked state attractors exist, but locked state attractors with 
more than three clusters are never possible 50 |. 



4. For C = 3, regions of parameter space exist where two frequency quasiperiodicity can 
occur with periodic attractors arbitrarily nearby in parameter space. 

5. For C > 4, (C — l)-frequency quasiperiodicity can occur, and periodicity and chaos 
occur for parameter values near those yielding (C — l)-frequency quasiperiodicity. 

This work is supported by the U.S. Army Research Office grant ^ W911NF-12-1-0101. 

Appendix A: Theoretical values of the critical coupling strength ^vith a uniformly 
distributed g{io) 

In this appendix we summarize the theoretical results of the critical coupling strength F^ 
when g{uj) is given by the uniform distribution 



^M = ^^(1-1^1), 



(Al) 



First, we determine F^ when there is no spread in /i(a), i.e., h{a) = 6{a — a), and 7 = 
for all oscillators. When a > 0, we have 



^ 1 TT 1 -1/1 

F-i = - + -tan-^ — 
4 2 \2a 



similarly, when a < 0, we have 



F, ^ = tan 



a\ 



(A2) 



(A3) 



For the cases when there is spread in h{a), we assume the same model h{a) = {26a) ^U{6a- 
\a — a\). By denoting a+ = a + 6a and a- = a — 6a, we have for a > 6a, 



^ 1 TT 1 

F7^ = - + 



4 4:6a 



a+ tan ^ 



2a, 



for a < 6a, 



r-' 



26a 



0:4. tan 



ftj 



«_ tan I — ^- \ -\ — In 
2a. 



-1 / M 1, 

a_ tan I — I + - m 

a- 




(A4) 



(A5) 
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and for lal < —5a, 



r 



-1 



1 1 + /2, where 



A = ?f^)+ ' 



4 \25a 

1 
'2fc 



4fc 



a+tan 



2a, 



-a_ tan 



— I — - In I a 
a. 



-ln|4ai + ll 



(A6a) 
(A6b) 

(A6c) 



Similar to the results with a Lorentzian g{uj)^ it can be readily shown that Eqs. (]A4p - (lA5p 
reduce to Eqs. flA2p and flASp in the limit (5a — )► 0. 

Next, we determine F^ when there is no spread in h{a), i.e., h{a) = 6{a — a) where a is 
constant, and the nonlinear frequency parameter 7 is a nonzero constant for all oscillators. 
For q; < 0, we know that 7 does not affect stability of the state (z) =0, so Fc is still given 
by Eq. f lASp . For a > 0, we have, by substituting s = iQ into the final expression after 
integration in Eq. f l25|) . and introducing ri± = Q±l + 0:7, that F^ and Q are to be given by 
the solution of the following pair of equations, 

V+iV- +4a^) 



c 



7 



= 7 



In 



IT — tan 



-1 fV± 
2a 



tan 



71 + tan 
-1 fV- 



-ifV± 
,2a 



2aJ\ 



1 



In 



tan-^ f ^ 
V2a 



ril{ril + Aa^) 



rj_{rj_ + 4a; 



M,2^ 



(A7a) 
(A7b) 



It can be easily checked from ( 1A7P that Fc reduces to ( lA2p when 7 — )• (Note fi — )■ in this 
limit) . 
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